MODULE fefmm_box_utils

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fefmm_init_box_utils,           &
             fefmm_shift_boxes_up_level,     &
             fefmm_shift_boxes_down_level,   &
             fefmm_get_boxes,                &
             fefmm_sh_pair_in_box,           &
             fefmm_box_centre,               &
             fefmm_box_map,                  &
             fefmm_get_neighbours,           &
             sefmm_get_filled_boxes,         &
             sefmm_AO_in_box
!             sefmm_get_parent_box,           &
!             sefmm_get_children

   REAL(REALK),   SAVE :: fefmm_coord_shift(3)
   REAL(REALK),   SAVE :: grain
   INTEGER(INTK), SAVE :: ndim(3)
   INTEGER(INTK), SAVE :: level

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_init_box_utils(grain0,sh_pairs)

      IMPLICIT NONE
      REAL(REALK),        INTENT(IN) :: grain0
      TYPE(fmm_sh_pairs), INTENT(IN) :: sh_pairs(:)

      REAL(REALK)   :: maxdim(3), mindim(3)
      INTEGER(INTK) :: i

      level = 1
      grain = grain0
      maxdim(:) = -1d99 
      mindim(:) = 1d99 
      DO i = 1, SIZE(sh_pairs)
         maxdim(:) = MAX(maxdim,(sh_pairs(i)%centre(:)+sh_pairs(i)%extent%hi))
         mindim(:) = MIN(mindim,(sh_pairs(i)%centre(:)-sh_pairs(i)%extent%hi))
      END DO
print *, 'maxdim', maxdim
print *, 'mindim', mindim
      fefmm_coord_shift = mindim
      ndim(:) = FLOOR( (maxdim(:)-mindim(:))/grain0 ) +1

   END SUBROUTINE

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_shift_boxes_down_level

      IMPLICIT NONE

      grain = half*grain
      ndim(:) = ndim(:)*2
      level = level + 1

   END SUBROUTINE fefmm_shift_boxes_down_level

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_shift_boxes_up_level

      IMPLICIT NONE

      IF (level == 1) RETURN
      grain = two*grain
      ndim(:) = ndim(:)/2
      level = level - 1

   END SUBROUTINE fefmm_shift_boxes_up_level

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_get_boxes(nbox)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(OUT) :: nbox

      nbox = ndim(1)*ndim(2)*ndim(3)
print *, 'ndim', ndim
print *, 'nbox', nbox

   END SUBROUTINE fefmm_get_boxes

!-------------------------------------------------------------------------------

   FUNCTION fefmm_box3c(box)

      INTEGER(INTK), INTENT(IN) :: box
      INTEGER(INTK) :: fefmm_box3c(3)

      fefmm_box3c(1) = MOD(box-1,ndim(1)) +1
      fefmm_box3c(2) = MOD(box-1,ndim(1)*ndim(2))/ndim(1) +1
      fefmm_box3c(3) = MOD(box-1,ndim(1)*ndim(2)*ndim(3))/(ndim(1)*ndim(2)) +1

   END FUNCTION fefmm_box3c

!-------------------------------------------------------------------------------

   FUNCTION fefmm_box1c(box)

      INTEGER(INTK), INTENT(IN) :: box(3)
      INTEGER(INTK) :: fefmm_box1c

      fefmm_box1c = box(1) + ndim(1)*(box(2)-1 + (box(3)-1)*ndim(2))

   END FUNCTION fefmm_box1c

!-------------------------------------------------------------------------------

   FUNCTION fefmm_box_map(box,x,y,z)

      INTEGER(INTK), INTENT(IN) :: box, x,y,z
      INTEGER(INTK) :: fefmm_box_map

      INTEGER(INTK) :: box3c(3)

      box3c = fefmm_box3c(box) + (/x,y,z/)

      fefmm_box_map = -1
      IF ( box3c(1) < 1       ) RETURN
      IF ( box3c(1) > ndim(1) ) RETURN
      IF ( box3c(2) < 1       ) RETURN
      IF ( box3c(2) > ndim(2) ) RETURN
      IF ( box3c(3) < 1       ) RETURN
      IF ( box3c(3) > ndim(3) ) RETURN
      fefmm_box_map = fefmm_box1c(box3c)
      !print *, box,';', fefmm_box3c(box),';', x,y,z,';', fefmm_box_map

   END FUNCTION fefmm_box_map

!-------------------------------------------------------------------------------

   FUNCTION fefmm_box_centre(box)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: box
      REAL(REALK) :: fefmm_box_centre(3)

      INTEGER(INTK) :: box3c(3)

      box3c = fefmm_box3c(box)
      fefmm_box_centre(:) = (box3c(:) - half)*grain + fefmm_coord_shift(:)

   END FUNCTION fefmm_box_centre

!-------------------------------------------------------------------------------

   FUNCTION fefmm_sh_pair_in_box(sh_pair,box)

      IMPLICIT NONE
      TYPE(fmm_sh_pairs), INTENT(IN) :: sh_pair
      REAL(REALK),        INTENT(IN) :: box(3)
      LOGICAL :: fefmm_sh_pair_in_box

      REAL(REALK) :: r12

! Select any box it overlaps with
      fefmm_sh_pair_in_box = .FALSE.
      r12 = half*grain + sh_pair%extent%hi
      IF ( ABS(sh_pair%centre(1) - box(1)) > r12 ) RETURN
      IF ( ABS(sh_pair%centre(2) - box(2)) > r12 ) RETURN
      IF ( ABS(sh_pair%centre(3) - box(3)) > r12 ) RETURN
      fefmm_sh_pair_in_box = .TRUE.

! Select just one box
!      fefmm_sh_pair_in_box = .FALSE.
!      r12 = half*grain
!      IF ( ABS(sh_pair%centre(1) - box(1)) > r12 ) RETURN
!      IF ( ABS(sh_pair%centre(2) - box(2)) > r12 ) RETURN
!      IF ( ABS(sh_pair%centre(3) - box(3)) > r12 ) RETURN
!      fefmm_sh_pair_in_box = .TRUE.

   END FUNCTION fefmm_sh_pair_in_box

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_get_neighbours(boxID,npts,hi,lo)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: boxID, npts
      INTEGER(INTK), INTENT(OUT) :: hi(3), lo(3)

      INTEGER(INTK) :: box3c(3)

!      lo(:) = npts+1
!      hi(:) = 2*npts
!return

      ! Identify only nearest neighbours not over the edge of grid
      box3c = fefmm_box3c(boxID)
      lo(:) = 1
      hi(:) = 3*npts
      IF (box3c(1) == 1)       lo(1) = npts+1
      IF (box3c(1) == ndim(1)) hi(1) = 2*npts
      IF (box3c(2) == 1)       lo(2) = npts+1
      IF (box3c(2) == ndim(2)) hi(2) = 2*npts
      IF (box3c(3) == 1)       lo(3) = npts+1
      IF (box3c(3) == ndim(3)) hi(3) = 2*npts
!      print *, 'lo =', lo
!      print *, 'hi =', hi

   END SUBROUTINE fefmm_get_neighbours

!-------------------------------------------------------------------------------

   FUNCTION sefmm_AO_in_box(grain,extent,AO_cntr,box_cntr)

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: grain, extent
      REAL(REALK), INTENT(IN) :: AO_cntr(3)
      REAL(REALK), INTENT(IN) :: box_cntr(3)
      LOGICAL :: sefmm_AO_in_box

      REAL(REALK) :: r12

      sefmm_AO_in_box = .FALSE.
      r12 = ROOT3*(half*grain) + extent

      IF ( ABS(AO_cntr(1) - box_cntr(1)) > r12 ) RETURN
      IF ( ABS(AO_cntr(2) - box_cntr(2)) > r12 ) RETURN
      IF ( ABS(AO_cntr(3) - box_cntr(3)) > r12 ) RETURN

      sefmm_AO_in_box = .TRUE.

   END FUNCTION sefmm_AO_in_box

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_get_filled_boxes(AOcntr,extent,box_list)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: AOcntr(3), extent
      INTEGER(INTK), POINTER, INTENT(OUT) :: box_list(:)

      INTEGER(INTK), TARGET, SAVE :: buffer(10000)
      INTEGER(INTK) :: box(3), tmp(3)
      INTEGER(INTK) :: i,j,k, ijk, width

      box(:) = INT(AOcntr(:)/grain) +1

      width = INT(extent/grain) +1

      ijk = 0
      DO i = -width,width
         tmp(1) = box(1) + i
         IF (tmp(1) < 1) CYCLE
         IF (tmp(1) > ndim(1)) CYCLE
         DO j = -width,width
            tmp(2) = box(2) + j
            IF (tmp(2) < 1) CYCLE
            IF (tmp(2) > ndim(2)) CYCLE
            DO k = -width,width
               tmp(3) = box(3) + k
               IF (tmp(3) < 1) CYCLE
               IF (tmp(3) > ndim(3)) CYCLE
               ijk = ijk +1
               IF (ijk > SIZE(buffer)) CALL fmm_quit('sefmm_get_filled_boxes')
               buffer(ijk) = fefmm_box1c(tmp)
            END DO
         END DO
      END DO

      box_list => buffer(1:ijk)

   END SUBROUTINE sefmm_get_filled_boxes

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_get_parent_box(box,parent)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: box
      INTEGER(INTK), INTENT(OUT) :: parent

      INTEGER(INTK) :: box3c(3), parent3c(3)
      INTEGER(INTK) :: nndim(3)

      box3c(:) = fefmm_box3c(box)
      parent3c(:) = INT((box3c(:)+1) / 2)

      CALL fefmm_shift_boxes_up_level
      parent = fefmm_box1c(parent3c)
      CALL fefmm_shift_boxes_down_level

   END SUBROUTINE sefmm_get_parent_box

!-------------------------------------------------------------------------------
!
!   SUBROUTINE sefmm_get_children(box,children)
!
!      IMPLICIT NONE
!      INTEGER(INTK), INTENT(IN)  :: box
!      INTEGER(INTK), INTENT(OUT) :: children(8)
!
!      INTEGER(INTK) :: box3c(3), tmp(3)
!      INTEGER(INTK) :: i,j,k, ijk
!
!      box3c(:) = fefmm_box3c(box)*2
!      CALL fefmm_shift_boxes_down_level
!
!      ijk = 0
!      DO i = -1, 0
!         tmp(1) = box3c(1) + i
!         DO j = -1, 0
!            tmp(2) = box3c(2) + j
!            DO k = -1, 0
!               ijk = ijk +1
!               tmp(3) = box3c(3) + k
!               children(ijk) = fefmm_box1c(tmp)
!            END DO
!         END DO
!      END DO
!
!      CALL fefmm_shift_boxes_up_level
!
!   END SUBROUTINE sefmm_get_children
!
!-------------------------------------------------------------------------------

END MODULE fefmm_box_utils
